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ABSTRACT 

A satellite galaxy or dark matter subhalo that passes through a stellar disk may 
excite coherent oscillations in the disk perpendicular to its plane. We determine the 
properties of these modes for various self-gravitating plane symmetric systems (Spitzer 
sheets) using the matrix method of Kalnajs. In particular, we find an infinite series 
of modes for the case of a barotropic fluid. In general, for a collisionless system, there 
is a double series of modes, which include normal modes and/or Landau-damped 
oscillations depending on the phase space distribution function of the stars. Even 
Landau-damped oscillations may decay slowly enough to persist for several hundred 
Myr. We discuss the implications of these results for the recently discovered vertical 
perturbations in the kinematics of solar neighborhood stars and for broader questions 
surrounding secular phenomena such as spiral structure in disk galaxies. 
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1 INTRODUCTION 

In the ACDM cosmological paradigm, galactic disks are em¬ 
bedded in extended halos of dark matter, which are pop¬ 
ulated by satellite galaxies, star streams and dark matter 
subhalos. Inevitably, some of this halo substructure will pass 
through the disk, heating and thickening the disk and also 
triggering the development of secular phenomena such as 
spiral structure, bars, and warps. 

In one of the earliest studies of disk heating, 
iToth fc OstrikeJ (Il992l) calculated the energy that a pass¬ 
ing satellite deposits in a stellar disk by assuming that the 
satellite gravitationally scatters indi vidual disk stars . Their 
analysis, which was based on the IChandrasekharl (Il943l ) 
dynamical friction formula, suggested that satellite infall 
could account for the thickness and velocity dispersion of 
galactic disks. On the other hand, the thinness and cold¬ 
ness of disks could be used to set constraints on the rate 
of satellite infall and hence the underlying cosmological 
model. Subsequent numerical experiments confirmed that 
satellites can in fact heat and thicken stellar disks, though 


triker (see, for example Ouinn et al. 

ll 19931): Walker et al.1 

lll996l): Huang & CarlbergI (Il997l): 

Benson et al.1 120041): 

iGauthier et al.1 (l2006l): IKazantzidis el 

al.1 120081)). 


portant effects. First, satellites are tidally disrupted by the 
gravitational field of the host galaxy, especially when their 
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orbits take them into the central disk-dominated region. Sec¬ 
ond, disks can respond coherently to the gravitational field 
of a satellite. Hence, much of the energy that is transferred 
to the disk takes the form of large-scale pert urbations such 
as warps in the outer disk llO uinn et"^ or a tilt of the 

disk plane (iHuang fc Carlberg|ll997l). Indeed satelli t es ca n 
change the morpholo gy of a galaxy. GjurthigreUal ll2006l) . 
iDubinski et al.l LoOSi ) , and iKazantzidis et ah 1 200g ) for ex¬ 
ample found that satellites and dark matter subhalos can 
provoke the formation of a ba r and/or spiral structure. More 
recently. IPurcell et ahl (1201 ll ) suggested that the Sagittarius 
dwarf spheroidal galaxy might be responsible for the Milky 
Way’s bar and spiral structure. 

The importance of coherent perturbations in disk - 
satellite interactions was stressed bv ISellwood et ahl lll998l) . 
They argued that a passing satellite transfers energy to 
the disk through the excitation of bending waves, which 
eventually decay through Landau damping. The upshot is 
that since the energy transfer from satellite to disk occurs 
through coherent perturbations disk heating is non-local. 


iToomrel (Il966l) derived the dispersion relation for bend¬ 
ing waves in a stellar disk of zero thickness and uniform sur¬ 
face density. In general, these waves propagate in the disk 
plane with a group velocity that depends on the properties 
of the disk and the wavelength of the perturbation. Though 
Toomre’s analysis showed that a perturbation with a suffi¬ 
ciently short wavelength is susceptible to a buckling insta¬ 
bility, he argued that this instability would be suppressed 
by the random motions of the sta rs in the vert ical direction. 
These results were confirmed by lAra kl lll985h who studied 
wavelike perturbations of a stellar disk with finite thickness. 
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In his equilibrium model, the surface density as well as the 

the Galactic midplane can 

result from satellite-disk interac- 

velocity dispersions in the vertical and horizontal directions 

tions dWidrow et al.ll2012l: 

Gomez et al.ll2013l:IWidrow et al.l 

are constant in space while the vertical structure is given bv 

l20l4l. Gomez et al.l (2013 

), for example, used N-body ex- 


lll95Ct) . He found that the system avoided the buckling insta¬ 
bility at all wavelengths provided the vertical velocity dis¬ 
persion wa s greater than 0.293 times the horizontal velocity 
dispersion. lArakil (1 19851') also considered breathing modes, 
which correspond to compression and expansion of the disk 
perpendicular to and symmetric about the disk midplane. 
In this case, the system can undergo a Jeans instability. 

Bending and breathing modes are the simplest modes in 
disklike systems, although they are no means only ones. In 
order to explore vertical perturbations in more d etail, one 
can ig n ore variations i n the horizonta l direc tion. lAntono-J 
lll97]I) : iKalnaid (Il973l ~): iFridman et al.l (Il984l ~). for example, 
considered strictly vertical perturbations of a homogeneous 
slab and found an inhnite doub le ser ies of normal mod es. On 
the other hand iMathuil (Il990ll and IWeinbei^ (Il99ll l inves¬ 
tigated vertical oscillations in a variant of the Spitzer sheet 
where the stellar energy distribution is truncated (the low¬ 
ered Spitzer sheet) and were able to identify only a handful 
of normal modes. 

Our aim is to carry out a comprehensive study of ver¬ 
tical oscillations in stellar a nd gaseous dis ks. I n particu¬ 
lar we extend the analysis of iMathurl ill99Cl[l and IWeinberd 
lll99lll to Landau-dampe d oscillations by using the matrix 
meth od of Kalnai^ (Il977l) and comple x anal ysis techniques 


from iLandaJ l 1946ll and iLynden-Belll lll962l 'l . We begin by 


calculating the normal modes of a gaseous disk, which serves 
as a warm-up to the the more complex collisionless case. We 
then discuss normal modes of the homogeneous slab. Finally 
we turn to the original (untruncated) and lowered Spitzer 
sheets. In both cases, we identify a double series of Landau- 
damped oscillations. We contend that t hese “modes”, com¬ 
bined with the normal modes found in IMathurl lll990ll and 
IWeinbe^ lll99lll . are analogous to the double series of modes 
of the homogeneous slab. Since we only consider vertical mo¬ 
tions we cannot directly address the question of how mode s 
propogate in the disk plane as was done in Arakil (Il985l l 


(see, also a brief discussion in IWeinberd (Il991 )) though we 


do gain a more complete understanding of vertical modes 
and Landau-damped oscillations. 

The primary motivation for this work comes from re¬ 
cent observations of vertical phase space structures in the 
kinematics of disk stars in the solar neighborhood of the 
Milky Way. These observations come free three surveys: 
the Sloan Extens ion for Galac t ic Un derstanding and Explo¬ 
ration (SEGUE; |Yanny__et_aL (l2()09|)) , the RAdial Velocity 


Experiment ('RAVE: lsteinmetz et al.l (l2006lH . and the LAM- 
OST Experime nt for Galac t ic Un derstanding and Explo¬ 
ration (LEGUE: iDeng et al.l (120121 ') L These surveys provide 
full six-dimensional phase space information for tens of thou¬ 
sands of stars within a few kiloparsecs of the Sun. Recently, 
several groups have de tected vertical bulk motions using 
data from these surveys (IWidrow et al.ll2012l : IWilliams et al.l 
I2OI3I : iGarlin et al.ll 20 f^ . which appear to take the form of 
compression a n d exp ansio n of the stellar disk . In ad dition, 
IWidrow et al.l (12012 1 and lYannv fc Gardn^ (l2013ll found 
evidence for wavelike North-South asymmetries in the num¬ 
ber counts of solar neighborhood stars. 

Velocity and number density perturbations normal to 


periments to show that a Sagittarius-like dwarf with a 
mass of 10^°'® — 10^*^ Mq could produce de nsity perturba¬ 
tions of the s ame amplitude as w as see n in IWidrow et al.l 


(l2012ll and lY annv fc GardneJ (l2013l l. Simulations by 
iFeldmann fc Spolvail (l2015l l showed that lower mass (10® — 
10® Mq) satellites produce subtle features in the bulk veloc¬ 
ity held of the disk that might be observed in the next gen- 
eration of astrometric surveys such as Gaia (IPerrvman et al.l 

I2OOIII . 

The bar and spiral structure of the Milky Way can also 
perturb the velocity distribution of stars in the disk and 
might therefore be responsible, at least in part, f o r the o bser¬ 
vation s des cribed above . In particular IPehnenl (l2000h . IFuxI 
(l200lll an dH ov 3(12013) showed that a resonant interaction 
between the Milky Way’s bar and stars in the solar neigh¬ 
borhood might be responsible for the Hercules stream, a 
group of co-moving stars whose b ulk velocity i s offset from 
that of the local standard of rest. iFaure et al.l (l20l4) used 
test-particle simulations to study the response of disk stars 
to a spiral potential perturbation and showed that spiral 
structure could generate vertical bulk motions in the disk 
akin to what has been observed in the SEGUE, RAVE, and 
LEGUE surveys. lDebattist3 (l2014f) carried out a fully self- 
consistent N-body simulation of a spiral galaxy and came 
to similar conclusions. Essentially, a spiral arm causes com¬ 
pression and expansion as it sweeps through the disk. 

Of course, since satellites have also been implicated in 
triggering the formation of bars and spiral structure, it may 
be difficult to disentangle perturbations that come directly 
from a disk-satellite interaction and perturbations due to 
structures in the disk, which themselv es were the result of 
passing satellites. IWidrow et al.l (l2014h foll owed the evolu¬ 
tion o f bending and breathing modes in the iGauthier et al.l 
(l2006ll simulation, where the disk was subjected to the con¬ 
tinual perturbations of a substructure-filled halo. Bending 
and breathing modes appear at early times as subhalos 
churn up the disk. The formation of a bar, which occurs 
at about 5 Gyr, is clearly triggered by substructure-disk in¬ 
teractions and it may well be that the breathing mode per¬ 
turbations are the mechanism by which this occurs. More¬ 
over, at late times the bar itself maintains both bending and 
breathing modes, especially in the inner parts of the galaxy. 

A second and more academic purpose for this work 
is to explore normal modes and Landau damping in self- 
consistent one-dimensional systems. In the usual textbook 
explanation of the Jeans instability one considers pertur¬ 
bations of a spatially homogeneous mass distribution. In 
the case of a fluid, one assumes that the sound speed 
is constant while for a collisionless system, one assumes 
a Maxwellian velocity distribution with constant velocity 
dispersion. In either case, one must confront the conun¬ 
drum that the unperturbed gravitational potential ipo is ill- 
defined. By symmetry, V'i/^o = 0 in a homogeneous system 
whereas Poisson’s equation implies S/'^'ipo = JttGpo. These 
two equations are inconsistent. The Jeans swindle, wherein 
one makes the ad hoc assumption that Poisson’s equation 
appl ies only to perturbed quant ities, provides a way forward 
(see iBinnev fc Tremaine! (l2008h for a more detailed discus¬ 
sion). The linearized equations are then solved by making 
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the ansatz that the perturbed density and potential vary 
harmonically in space and time. In both the fluid and col¬ 
lisionless cases a perturbation whose wavelength is longer 
than the Jeans length grows exponentially. A perturbation 
in a fluid whose wavelength is less than the Jeans length 
oscillates as sound waves. On the other hand, a short wave¬ 
length perturbation in a collisionless system undergoes Lan¬ 
dau damping and rapidly decays. In the models considered 
here, the density is concentrated in the midplane and the 
Poisson equation can be solved without resorting to the 
Jeans swindle. 

The outline of the paper is as follows: In Section 2 we 
find the normal modes for an isothermal plane-symmetric 
fluid. In Section 3, we derive normal and Landau-damped 
oscillations for the three examples of collisionless systems 
mentioned above. In Section 4, we present results from a sim¬ 
ulations of a one-dimensional system that exhibits damping 
in a perturbed Spitzer sheet. We conclude in Section 5 with a 
summary and discussion of our results and some thoughts on 
further directions for this line of research. Three appendices 
provide mathematical details for some of our calculations. 


has an equation of state p = p{p) and constant sound speed 
Vs = {dp/dpf^^. 

We write p, v, p, and ijj as the sum of an equilibrium 
solution and a linear perturbation, e.g., 

p^ Po{z) + pi{z,t) . (4) 

For the equilibrium solution, vo = 0 and the continuity equa¬ 
tion is satisfied automatically while the Euler and Poisson 
equations are solved by the following density-potential pair: 

V’o (z) = 2vl In (cosh {z/ zq)) (5) 

and 

po (z) = Pc sech^ {z/zq) (6) 

where zq = Vs! {2 'kGpcV^^ and P e is the density in the mid¬ 
plane (ISDitzeiill942l : [^mmlll95(]ll . 

In what follows, we use a system of units in which 20 = 
Vs = G = 1 and pc = l/27r. The linearized equations are 


2 LINEAR PERTURBATIONS OF AN 

ISOTHERMAL ONE-DIMENSIONAL FLUID 


dpi d . ^ „ 

^ ^ (po^i) = 0 


( 7 ) 


In this section and the ones that follow, we consider the ver¬ 
tical perturbations of plane symmetric systems. The idea of 
treating a galactic disk as a pl ane symmet ric system can be 
traced to the seminal paper bv lOort lll932l L wherein the mo¬ 
tions of stars perpendicular to the Galactic midplane were 
used to estimate the verti cal force and mass distribution 
in the solar neighborhood. ISoitzeil (Il942l l derived an equi¬ 
librium model for a self-gravitating, plane-symmetric sys¬ 
tem of stars under the assumption that the stellar velocity 
dispersion is constant with height above the midplane. His 
derivation is based on the Jeans equations (i.e., moments of 
the collisionless Boltzmann equation) and is therefore akin 
to the derivati on of the equi librium fluid model considered 
in this section. ICamml Jl950ll solved for the equilibrium dis¬ 
tribution function, which provides the starting point for our 
analysis in Section 3. In what follows, we refer to these mod¬ 
els collectively as Spitzer sheets. 

We consider linear perturbations in a simple self- 
gravitating barotropic fluid. For an analysis of perturba¬ 
tions in a m ultiphase, magnetized model of the interstellar 
medium see lWalters fc ll200ll L A fluid with density p, 
velocity u, pressure p, and gravitational potential tp obeys 
the continuity, Euler, and Poisson equations. For a plane 
symmetric system, these equations become 


dvi dpi dipo dipi 

dt dz dz dz 


(8) 


d'^ipi 

dz^ 


inpi . 


( 9 ) 


To search for modes, we assume that the first-order quan¬ 
tities are proportional to exp{—iujt) (e.g., pi{z,t) = 
exp (—ituf) pi( 2 )) and combine the continuity and Euler 
equations to arrive at a single equation for pi: 


2 - 


VJ pi = - 


d^p'i 

dz^ 


dpi dipo 
dz dz 


dz^ 


dpo dipi d^ipi 

dz dz dz'^ 


( 10 ) 

( 11 ) 


For a homogeneous system, the fourth term on the 
right-hand side is zero while the second and third 
terms are ignored by im plementing the Jeans swindle 
dBinnev fc Tremaind l2008l ~l . We then assume a sinusoidal 
spatial dependence for the perturbation (pi, ipi oc exp (ikz)) 
and use the first-order Poisson equation to arrive at the dis¬ 
persion relation 


dv dv _ I dp dip , , 

dt ^ dz p dz dz 

g = m 

where we choose our coordinate system so that the 2 -axis 
is normal to the symmetry plane. We assume that the fluid 


uF = Vs (D -kj) =D -2 (12) 

where kj = [AnGpc/vf) = 2^^^ is the Jeans wavenumber. 
For k < kj (long wavelength perturbations) oj is imaginary 
signaling an instability with growth rate a = (2 — 

On the other hand, for k > kj perturbations oscillate with 
frequency cu = (fc^ — 2)^^^. 

To study linear perturbations of the Spitzer sheet 
we use the Kalnajs matrix method and write pi and ipi 
in terms of a biorthono rmal basis llKalnaisI Il97ll . Il977l : 
iBinnev fc Tremainell2008ll : 
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^1 G) = cj v-i ( 2 ) pi ( 2 ) = 

(13) 

where 



(14) 

and 


POO 

/ dz'ip* {z)pk{z) = -Sjk ■ 

J — OO 

(15) 

We multiply Eq.[Tn]by — and integrate with respect to 2 
from — OO to 00 to obtain a matrix equation of the form 

LO Cj -— ^ ^ AdjkCk . 
fc 

(16) 

Evidently, the eigenvalues of M are the squares of the mode 
frequencies. 

For the calculation at hand we use the basis introduced 
inlArakil (1198,4): 

i’j G) = EjPj{u) 

(17) 

and 


Pi G) = (1 - u^) p,{u) 

(18) 


where u = tanh( 2 ), Pj are Legendre polynomials, and Nj = 
{2n {2j + 1) /j {j + are normalization constants. The 

matrix elements Mjk can be calculated analytically (see 
Appendix A). Note that Mjk is nonzero if and only if 
j = k — 2, k, OT k + 2. The matrix therefore separates into 
two independent matrices, one where j and k are both even 
and the other where j and k are both odd. The physical im¬ 
plication is that modes have definite parity, a property that 
follows from the symmetry of the equilibrium system under 
the transformation 2 —>■ — 2 . 

The eigenvalues and eigenvectors for these sparse ma- 
trices are found usin g the software package LAPACK 
ll Anderson et al.lll992l l. The frequencies uin and eigenfunc¬ 
tions for the lowest eight modes are given in Table 1 and Fig¬ 
ure [T] The eigenfunction label n corresponds to the number 
of nodes in the density profile. Note that the n = 1 mode, 
which has zero frequency, corresponds to a shift in the sys¬ 
tem as a whole. 

3 COLLISIONLESS SYSTEMS 
3.1 Formalism 

The dynamics of a collisionless, plane symmetric system is 
described by the collisionless Boltzmann and Poisson equa- 
tions in one d imen sion. We follow the formalism found in 
iMathurl (Il990l ~l and IWeinberd (Il99ll ') in which the collision¬ 
less Boltzmann eq uation is written in terms of angle-action 
variables (see, also lBinnev fc Tremaine! ll2008l ) ). For an alter - 
nate approach based on Jeans equations, see iLoui i (Il992ll . 
A particle in a time-independent potential, ipoG) executes 
periodic motion with constant energy E = /2 + %l)o {z), 

period T{E), and maximum excursion from the midplane 
2max where tpo (zmax) = E. We can therefore introduce the 



012345012345 
Z Z 


Figure 1. Normal modes of the fluid Spitzer sheet. Panels on the 
left show the potential (top) and density (bottom) for the odd 
parity modes with n = 1 (red curve), n = 3 (green), n = 5 (blue) 
and n = 7 (magenta). Panels on the right show the even parity 
modes n = 2 (red curve), n = 4 (green), n = 6 (blue) and n = 8 
(magenta). 


angle-action variables (6^, E) where 0 is defined so that dt = 
T{E)d9/2TT and 6 — 0 corresponds to {z, v) = (—2max, 0). 

As before, we write the density and potential, along 
with the distribution function, as the sum of an equilibrium 
solution and a linear perturbation. For example 

/ {E, 9, t) = /o {E) + /i (E, 6, t) . (19) 

The linearized collisionless Boltzmann and Poisson equa¬ 
tions are then 

dfi Stt dfi _ 2n dfp dtpi ^ 
dt T{E) 89 T{E) dE 89 '' ’ 

and 

= 47vG j dv fi . (21) 

Note that the Newtonian potential that appears in the col¬ 
lisionless Boltzmann equation includes contributions from 
both the system and any external perturbations while fi 
refers only to the system. 

Since particle orbits are periodic in 9 we can expand 
each of the first-order quantities in a Fourier series, e.g., 

OO 

h{E,9,t)= UiE,t)E^\ (22) 

n= —OO 

Eq. [20]is satisfied for each Fourier component: 
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odd modes 


even modes 


n 


n 

^71 

1 

0 

2 

1.129 

3 

1.522 

4 

2.180 

5 

3.138 

6 

4.373 

7 

5.976 

8 

7.865 


Table 1. Frequencies, in units of where Vs = zq = 1 for the first four even and first four odd modes. 


dfn , 2'Kin j. 2-Kin dfo ^ 
dt T{E) dE'^" ^ 

The temporal Fourier transform of f„ is 

roo 

U{E,w)= dtf„{E,t)V^^ 

Jo 


(23) 

(24) 


and similarly for ■(/)„. In writing /„, we have assumed that 
/„ = 0 for t < 0 and that there exists a real number 7 > 0 
such that f dtexp(—yt)f„(t) converges. The latter condi- 
tion insures the exis t ence o f the inverse Fourier transform 
dBinnev fc Tremain3 (120081 ') - see Eq.[35] below. We multi¬ 
ply Eq.[^ by exp (icut) and integrate over t. The first term 
is handled by an integration by parts and we obtain an al¬ 
gebraic equation for /„, which can be written as 


_ dfo nQ{E}ip„ (E, lj) 
dE nfl{E) — uj 


(25) 


where fI(E) = 2k/T{E) is the orbital frequency for a parti¬ 
cle with energy E. 

As in the fluid case, we decompose the density and po¬ 
tential in terms of a biorthonormal basis. Let tpi be the 
linear perturbation to the potential due to the system itself 
and i/if be the linear external potential (if one exists). For 
the system, we have 


00 

i’l 0, w) = ^ Cj (uj) ipj (z) (26) 

j=i 

and 

/ CO 

dvfi=Y^ Cj (w) pj (z) (27) 

where the potential-density pair {tpj, pj) satisfy Eas. ll3lfT^ 
Likewise, for the external potential, we have 


f’l 0, ^) = Si V) ■ (28) 

We then write ipj as a Fourier series 


00 

Yf ^UE)V^^ (29) 

n= —00 


We now use Eq.[25]and the expansions defined above to write 
the Poisson equation as follows: 


1 

4:kG dz'^ 


CO 


y ] Pi 



dfo n^-ipn 

-—e 

dE nil — LO 


(31) 


where = J2j Vi + dj) 4>jn. 

We multiply this equation by —ifj (z) and integrate over 
2 to arrive at a matrix equation the expansion coefficients 
Cj- 


Cj (uj) = Mjk (uj) (cfc (w) -I- dk (uj)) (32) 


where 


Mjk (uj) = -4k 



iTj^dfo ^f^jn'^kn 
-tj 2 /n 2 


(33) 


(^Mathul^ll990^ : IWeinberdll991^ . Note that in deriving Eg . 1331 
we have combined positive and negative values of n and 
also omitted the n = 0 term as its contribution to the 
sum is zero. The matrix M, which is sometimes referr ed to 
as the polarization matrix (iBinnev fc Tremaine! 1200^ . de¬ 
scribes the response of the system to the total perturbing 
potential (system plus external). We can rewrite Eg . 1321 as 
an explicit equation for the Cj-. 


Cj {^) = Y -^ 1 '“ 

k 

where (I-M)"^M describes the response of the sys¬ 
tem to an external potential. 

The time-dependent response of the system is deter¬ 
mined by calculating the inverse Fourier transform of Cj: 

1 ri^+oo 

Cj (f) = ^ 3 - / du Cj (uj) (35) 

= —i 'Y^ Cj^r (kJr) (36) 

r 

where uir are the poles of the function Cj (uj). These poles 
occur at the poles of the response matrix R or equivalently, 
where either det(M) is singular or det (I — M) = 0. The 
second of these conditions is equivalent to the eigenvalue 
equation 


where 

1 

^r.j{E) = — dOfjj {z {E, e)) e-™'’ . (30) 


cj (t*^) = Y -^ 1 * 

k 

Our goal is therefore to find frequencies tOm such that one 
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of the eigenvalues of M is equal to unity. The corresponding 
functions ■i/'i G, ^m) and p\ [z, uJm) are then the potential- 
density pairs for either normal modes, when Im {tOm) = 0 or 
Landau-damped oscillations, when Im {uJm) < 0. 


3.2 Homogeneous slab of stars 

We first deter mine the norma l modes of the spatially homo¬ 
geneo us slab (lAntonovlll97ll : lKalnais|ll973l : iFridman et al.l 
I1984 1. where the distribution function for the equilibrium 
system is given by 


,38) 

I 0 otherwise . 

The density is equal to the constant pc = for 

l-^l < /2'kGpc)^^^ and zero otherwise. All particles ex¬ 

ecute simple harmonic motion about the midplane with fre¬ 
quency flc = (dTrGpc)^^^. We choose units akin to those used 
in the previous section, namely a — G — 1 and pc = 1/2 tt. 
The edge of t he slab is then at z = ±1 and flc = 2^^^. 

Following lKalnai3 (Il973l i we use the biorthonormal ba¬ 
sis to describe the density and potential inside the slab: 

tpjiz) = Nj iPj+i{z) - Pj-i{z)) (39) 


and 


k 


o 




j 


Figure 2. Frequencies for normal modes of th e homogeneous 
stellar slab. The dots give the frequencies from iKalnai J <ll973ll 
for different values of the eigenvalues j and k. Red dots are for 
even j and correspond to modes where the density and potential 
have odd parity under 2 —y —z. Blue does are for even parity (odd 
j) modes. Modes in the same column have the same density and 
potential but different phase space distribution functions. 


pAz) = 


Mjk (uj) 


j (42) 


where Nj = {2ti/ (2j + 1))^^^. In fact, (ipj, pj) are the nor¬ 
mal modes of the system. We see that the density and 
potential have even parity when j is odd and vice versa. 
Further more, when j is odd, there are {j + 1) /2 distinct 
modes jKalnai d Il973l l . We label the frequencies for these 
modes LJj^k where k = 1, 3 ... (ji + 1)/2. Likewise, when 
j is even, there are {j -|- 2) /2 modes. For example, there 
is a single mode with j = 1: '01 = {'i N\/2) (z^ — l) . 
Pi = SNi/^tv , and oji.i = 6 ^'^^ = 1.73flc il Antonov! 1 197ll : 
iKalnaisI 1 19731 : [Fridman et al.lll984l L The distribution func¬ 
tion for this mode is 


SNiE dfo cos 2^ 

2 dE2-Lu^/4: 


(41) 


(see Appendix B). The frequencies of the j ^ 5 mode s are de¬ 
picted in Figure[2] numerical values can be found in lKalnal^ 
lll973ll . 


Self-consistency requires pi = J dv fi. Formally, this 
integral is inhnite because of the strong divergence of dfo/dE 
as J5 —>■ 1. To handle this divergence we consider a family of 
functions f\{E) that are continuous, that vanish as —>■ oo 
and that approach Eq.[3Hliii the limit A —>■ oo. We can then 
perform an integrati on by parts an d after wards let A —^ 
oo (see for example, [Fridman et al.l lll984h i. As shown in 
Appendix B, this procedure leads to the required relation 
between pi and fi. A similar trick can be used to evaluate 
the matrix elements in Eq. [33] where we can write 


The function fo has an integrable singularity at E = 1 and 
the matrix elements can be calculated numerically. In do- 
ing so, we hn d that M is indeed diagonal, a result that 
iKalnaisI lll973ll showed using various relations between Leg¬ 
endre polynomials and hypergeometric functions. 

In Figure [3] we show the density and potential pairs for 
the even parity modes with j = 1, 3, 5. The case j = 1 is a 
breathing mode in which the density within the slab depends 
on t but not z, while the boundaries of the slab move in 
and out accordingly. The distribution functions for the six 
distinct modes with these values of j are shown in Figure ID 
Modes in the same row have frequencies that cluster about 
{j -I- 1) Qc, as in Figure [2] For example, the three modes 
along the bottom row have frequencies from left to right of 
Lo/^c = 1.73, 2.07, 2.01. Particles make roughly half an orbit 
in the time it takes the system to execute a single mode 
oscillation of this type. The frequencies for the two modes 
in the middle row are uj/Q.c = 3.83, 4.06 so that particles 
make roughly one quarter of an orbit during one oscillation 
period. Modes in the same column have the same density 
profile and potential since they have the same eigenvalue j. 


3.3 Linear perturbations of the stellar Spitzer 
sheet 

In a collisionless system, oscillations that are in resonance 
with any of the particles in the distribution will decay via 
Landau damping whereby coherent energy in the oscillations 
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Figure 3. Even parity normal modes of the homogeneous slab. 
The linear potential (top panel) and density (bottom panel) are 
shown as a function of z for j = 1 (red/solid), j = 3 (blue/dotted) 
and j = 5 (green/dashed). 




Figure 4. Distribution functions for even parity modes of the 
homogeneous slab with j = 1, 3, 5. The arrangement of the panels 
is the same as the arrangement of the blue dots in Figure[2l Thus, 
the three modes along the bottom row have, from left to right, 
j = 1, 3, and 5 and uj ~ 1.73, 2.07, and 2.01. The three modes 
along the right column all have j = 5 and have a density profile 
and potential given by the green/dashed curves in Figure O 


is transferred to energy in random particle motion. The ho¬ 
mogeneous slab is a special example where all particles have 
the same orbital frequency fie- Since the mode frequencies 
are all slightly off resonance there is no Landau damping, at 
least not formally. 

In this section we consider the stellar Spitzer sheet 
where the distribution function is 

fo{E) = . (43) 


Figure 5. Range of allowed frequencies for orbits of stars and 
their ha rmonics in the low ered Sp i tzer s heet (LSS), the Spitzer 
sheet of ISpitz^ <ll942l) and ICamrd lll950t) (SS) and the homoge¬ 
neous plane (HS). The left-most horizontal bar in the row marked 
LSS shows the range of allowed frequencies for orbits in the low¬ 
ered Spitzer sheet with W = 2. The blue bar immediately to the 
right of this shows the frequencies for the second harmonics of 
these orbits, and so o n. The vertical arro ws mark the pos itions 
of the modes found by iMathuil lll99(ll and IWeinber^ lll99ll 'l . The 
frequency ranges for the Spitzer sheet are displayed in the same 
way while for the homogeneous plane all particles orbit at the 
same frequency. Odd parity modes and frequencies are shown in 
red while even parity modes and frequencies are shown in blue. 


The density and potential have the same form as in the 
fluid case (Eqs.[S]and|3 but with Vs replaced by a and /c = 
pej . As before, we choose nnits in which G = zo = 

(J = 1. In these units, the orbital frequencies range from 
Q{E) Qc — 2^/^ for i? ^ 0 to Q{E) 0 for E —>■ 
oo. Thus, coherent oscillations with frequency to will be in 
resonance with particles whose orbital frequencies are fin = 
u}/n where n = Timin, u-min -I- 1, . .. and rimin = Int(f2c/i^). 
This point is illustrated in Figure O The line segments in 
the middle portion of the diagram show the range in orbital 
frequencies and higher harmonics for the Spitzer sheet. The 
dots along the top portion of the diagram show the same 
for the homogeneous slab. The modes of the homogeneous 
slab are all off-resonance, that is, reside in the gaps between 
the dots. No such gaps exist in the case of the Spitzer sheet. 
The implication is that all coherent oscillations of the Spitzer 
sheet are damped. 

We also consider the lowered Spitzer sheet, which pro¬ 
vides a link between the the homogeneous slab and the 
Spitzer sheet. This model the one-dimensional analog of the 
lowered isothermal sphere or King model and has a distri¬ 
bution function that is given by 


ME) = 


0 


G < E <W 

otherwise . 
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With a = 1 and 


fw = cr^ (^(2Tvf^^ erf (yo) - 2®/^7ryoe (45) 

where yo = the density in the midplane is pc = 

l/27r and the maximum orbital frequency is flc = 2^^^ (See 
Appendix C2 for details). For finite W there is a lower bound 
on the particle frequencies: flmin = 0,{W). 

The horizontal line segments in the lower portion of 
Figure[5]show the range in frequencies and higher harmonics 
for the case W = 2. We se e that there a re n o resonant 
partic les for 0 < a; < 0.72 flc. iMathuJ (Il990[ l and IWeinbei^ 
lll99ll l refer to this range in frequency as the principal gap. 
Likewise, there are no particles with < lo < 1.45 flc (the 
hrst gap). 

As discussed above, our aim is to find values of to for 
which one of t he ei genvalues o f the m atrix M is equal to one. 
iMathurl lll990h and Weinberg (Il99ll ~) restrict their search for 
modes to real values lu that reside in the frequency gaps. In 
doing so, they avoid any potential singularities in the en¬ 
ergy integrals necessary for calculating the matrix elements. 
The vertical arrows in Figure [^indicate the posi t ions o f the 
modes found by the method outlined in lMathu3 (ll99Clfl and 
IWeinbe^ lll99 ih . Note that apart from the cu = 0 mode, the 
frequencies are slightly less than nflmin. 

Next we consider complex values of to. Recall that the 
Fourier transform /„ and hence the matrix elements Mjk are 
defined for Im(tj) > 0. In this regime, the energy integrals 
in Eq.[33]can be calculated by a straightforward numerical 
integration since the singularity in the integrand occurs in 
the lower half of the complex E plane and is avoided as one 
integrates along the real E axis. 

To explore solutions with Im(cu) ^ 0 we analyt¬ 
ically continue M into this region of the complex oj 
plane. To see how this works we dehne Mjk {oJr, ^i', A) = 
Mjk Gr + iXoJi) to be a function of the real variable A with 
LOR and Loi taken to be positive constants. The analytic con¬ 
tinuation of Mjk will have the property that Mjk is a con¬ 
tinuous function of A. When A changes from positive to neg¬ 
ative values, the singularity in the energy integral crosses 
from the lower half of the complex energy plane to the up¬ 
per half. Thus, if Mjk is to be continous when this hap¬ 
pens, the integration contour for the energy integral must 
deformed so that the It always remains above the singularity. 
This deformed contour is analogous to the so-called Landau 
contour used to carry out the velocity-space integral in the 
usual derivation of Landau damping fo r a homogeneous sys- 
tem (lLandaulll946l : lLvnden-Belllll96^ : iBinnev fc Tremainel 
l2008h . We then have two contributions to the energy inte¬ 
grals in Eg. 1331 one from the principal part of the integral, 
and the other from the residue due to the singularity at 
Q,{E) = Lo/n. That is, each term of the sum in Eg . 1331 be- 
comes 


— AttV 


dE 


^'^j‘n‘4^kn dfo 2^ 
-uj^/n^dE^ "" * 




~^dE 


Q=tj / n 

(46) 


The second term requires us to treat the quantities inside the 
parenthesis as complex analytic functions of complex energy. 
Details of how this is done can be found in Appendix C. 
We evaluate M as a function of complex to for 0 < 
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Figure 6. Heat map showing the position of (Landau) modes in 
the complex frequency plane for the lowered Spitzer sheet with 
W = 2 (top panel) and W = 5 (middle panel) and for the isother¬ 
mal plane (bottom panel). Colors indicate log-base 10 of the eigen¬ 
value of the matrix Mij — 5ij that comes closest to zero. 


Re(a;) < 8 and —0.6 < Im(a;) < 0. Eor e ach u we deter¬ 
mine the eigenvalues of M using LAPACK llAnderson et al.l 
I1992II and Hnd the eigenvalue closest to unity. A* Eigure [6] 
shows heat maps of 77 = log(|A* — 1|) as a function of lo 
for a lowered Spitzer sheet with W = 2 and W = 5 and for 
the infinity (i.e., W 00) Spitzer sheet. Places where H 
is a large negative number indicate mode frequencies. The 
matrix M is truncated at j =8, which explains why only 7 
or 8 candidate modes are found. In the W = 2 case, there 
is a single undamped breathing mode with lo ~ 1.4IIc as 
well as seven damped oscillations. For the latter, the ratio 
of the exponential decay time to the oscillation period is 
r/T = Re(a))/27rlm(a;) ~ 2.5 — 8. 

Figure |6] focuses on modes with frequencies clustered 
around Re (w) = 2Qc but with very different spatial struc¬ 
tures. Thus, they would correspond to the left-most column 
in in Figure 4 or the column of modes at a; ~ 2 in Fig¬ 
ure 2. In Figure [7] we zoom out for a wider view of the 
imaginary-Lu half-plane to show positions of modes clustered 
about Re (w) /Dc = 2, 4, 6. 


4 SIMULATIONS 

In this section, we present results from a simple one¬ 
dimensional simulation that illustrates the behaviour pre¬ 
dicted by our linear theory calculations. Our N-body sys¬ 
tem comprises self-gravitating, collisionless inhnite sheets. 
We use a particle-mesh scheme in which the density is cal¬ 
culated on a one-dimensional grid and the force on a sheet 
at position z is determined from the integral 

F{z) = 2'kG f dz p{z )sgn — z) . (47) 
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Figure 7. Positions of solutions to the dispersion equation for 
the W = 2 lowered Spitzer plane. The cluster of modes at 
Re(a;)/Oc ci 1.7 and Im(a;)/Oc ~ —0.07 is the same as the modes 
depicted in the top panel of Figure |6] 

Similar simul a tions were presented in Wember'3 lll99lll : 
IWidrow et al.l l|2012h and IWidrow et al.l ( 2014l 'l. 

We use initial conditions that correspond to a simple 
breathing mode by first setting up an equlibrium distribu¬ 
tion and then perturbing the positions and velocities accord¬ 
ing to the relations z = XzZe and v = XvVe where (^e, Ve) are 
the phase space coordinates of a particle in the unperturbed 
system and Ae,u are constants. With our choice of units, 
the total kinetic and potential energies of the unperturbed 
system are T = l/2n and V = 1/tt respectively. With this 
in mind, we choose = (3 — 2A2)^^^. The total energy of 
the perturbed and unperturbed systems are therefore the 
same while the initial virial ratio for the perturbed system 
is R = 2T/V = i3-2Xz)/Xz. 

In Figure[8]we show the time evolution of various system 
properties for a simulation with 2QQK particles and Xz = 
0.9. For example, in the upper left panel, we show the virial 
ratio R. During the initial phase from f = 0 to f ~ 12 
the amplitude of the oscillations in R damps from an initial 
value of 0.3 to ~ 0.02. The time-dependence of R during this 
phase can be described by a model comprised of the sum of 
damped exponentials: 

n 

R = Ro+ ^ Ci cos {uJit - . (48) 

i=l 

The fit for the run shown in Figure [H] has two terms with 
(tUi, a?i) ~ (1.5, 0.25) , (2.0, 0.24). Note that we only fit 
Eq.|48]for t < 12. At later times, perturbations are driven by 
the Poisson noise of the simulation (see below). The simula¬ 
tion results are in good qualitative agreement with our ana¬ 
lytic results. In particular, the lower panel of Figure [6] shows 


Figure 8. Time evolution of a perturbed plane symmetric col¬ 
lisionless system. The upper left panel shows the virial ratio 
R = 2TjV as a function of time (black solid curve) as well as 
a two-term fit based on Eg. 1481 (red dashed curved). The lower 
left panel shows the first four even coefficients j (j + 1) Cj that 
appear in the expansion for the density. Colors are red, blue, ma¬ 
genta, and green for j = 2, 4, 6 and 8 respectively. The upper 
right panel shows the time evolution of the width of the region 
that contains the inner 50% of the mass. The lower right panel 
shows the same for the inner 90% of the mass. 


that the typical modes of the Spitzer sheet have an oscilla¬ 
tion frequency of Re(a;) ~ 2 and an exponential damping 
constant of Im(a;) ~ 0.2. However, because the oscillations 
damp after only a few cycles we are unable to resolve the 
detailed mode structure anticipated in Figure [O] 

A key feature of vertical oscilations is that they damp 
most rapidly in the inner parts of the distribution. This is il¬ 
lustrated in the lower left panel of Figure [8] where we see that 
the damping rate decreases with increasing j. Furthermore, 
if we compare the top and bottom panels on the right, we 
see that the oscillations of the region containing 50% of the 
mass damp more rapidly th an the oscilatio n s of t he region 
containing 90% of the mass. IWidrow et al.l (l2014l l also no¬ 
ticed that the high energy part of the distribution function 
was more susceptible to vertical oscillations that were pro¬ 
voked by a passing satellite. We can therefore predict that 
bulk vertical motions will be strongest among stars in the 
high (vertical) energy tail of the phase space distribution. 

At late times, the initial perturbation has damped away 
but oscilations, seeded by Poisson noise due to the finite 
number of particles, persist. In Figure [9] we show late time 
behaviour of a low-resolution (lOK particle) simulation that 
was evolved until t ~ 1300. We show a portion of the 
time evolution of the virial ratio (upper panel) and the 
power spectrum of the time-domain Fourier transform (lower 
panel). We see that there are continous oscillations with ran¬ 
domly changing phase and amplitude but with a character¬ 
istic frequencies of cj ~ 2 and cj ~ 4. 
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Figure 9. Results from a low resolution simulation with lOK 
particles that was run until t = 1300. The upper panel shows a 
small section of the time evolution of the virial ratio while the 
lower panel shows the power spectrum. 


5 DISCUSSION AND CONCLUSIONS 


A passing satellite can irrevocably change the properties of a 
galactic disk in two fundamental ways. First, orbital energy 
from the satellite can be transferred to disk stars thereby 
heating and thickening the disk. Second, the satellite can 
provoke the formation of secular phenomena such as a bar, 
warp, or spiral structure. In either case, the initial response 
of the disk to the satellite can involve coherent oscillations. 
The subsequent evolution of these motions is governed by 
a mix of processes that include the restoring force of the 
disk’s own self-gravity, Landau damping, differential rota¬ 
tion, and swing amplification of spiral waves. In this paper, 
we have focused on the first two of these effects by restrict¬ 
ing our attention to plane symmetric systems. Indeed, the 
structure of the oscillatory behaviour for a self-gravitating, 
plane-symmetric system is already quite complicated, es¬ 
pecially in the case of a collisionless (i.e., stellar) system. 
There, the modes form a double series defined by two eigen¬ 
values, one that determines the mode’s frequency and the 
other that determines its spatial structure 

Linear perturbations of both gaseous and stellar plane 
symmetric systems divide neatly according to their parity 
with respect to the Galactic midplane. Th e North-South 
asymm etri es in the number counts found in IWidrow et al.l 
(l2012ll and lYannv fc Gardner! ll2013h by construction picked 
out odd parity density modes. Bending modes, which cor¬ 
respond local displacements of the disk from the midplane, 
can also be thought of as odd parity density perturbations. 
The midplane displaceme nts seen in the interstellar medium 
llNakanishi fc Sofuell2006n are an example of this. 

By contrast modes where the density perturbation has 
even parity have bulk velocity fields that are odd in parity. 
The simplest example is the breathin g mode perturbation , 
which may be present in the SEGUE (I Widrow etjah 201^ 
RAVE (IWilliams et Id] 1201311 . and LAMOST (iGarlin et al.l 


I 2 OI 3 I) data sets. Unfortunately, a clear picture of the ve¬ 
locity field throughout the 1-2 kiloparsec neighbourhood of 
the Sun is lacking in large part because of the complicated 
and in complete footprints of the surveys (see iGarlin et al.l 
ll2013ll for a detailed discussion). A combined analysis of the 
three data sets might improve this situ ation as will satellite 
data f rom Gaia dPerrvman et al.ll200ll l. As well. lBovv et al.l 
ll2014h proposed an alternative way to characterize the ve¬ 
locity field of the Galactic disk. The idea is to calculate the 
power spectrum of the velocity field after subtracting an ax- 
isymmetric model that accounts for the rotation of the disk. 
They find that the largest contribution to the power spec¬ 
trum on large scales comes from the Sun’s motion relative 
to the local standard of rest. In addition, they find a broad 
peak in the power spectrum with 0.2kpc“^ < k < 0.9kpc“^ 
and the associated motions might be associated with the 
time-dependent gravitational potential of the bar. While 
their analysis uses only heliocentric line-of-sight velocities 
and focuses on the two-dimensional in-plane velocity field, 
the technique could easily be extended to include proper 
mot ions and velocities perpendicular to the Galactic plane. 

iRobin et al.l ll2003h estimate the density of stars and 
dark matter in the midplane of the Galaxy at the position 
of the Sun is pc — 0.055 Mqpc“®. The period of vertical 
oscillations for a star near the midplane is therefore Tc — 
llOMyr or approximately one half the period for a star at 
this radius to orbit about the center of the Galaxy. In Section 
3, we found that in the solar neighbourhood, the ratio of the 
exponential decay time to the vertical oscillation period for 
pure vertical modes typically fell in the range r/T ~ 2.5 — 8. 
Thus, these modes might be expected to persist for 1—4 
orbital periods or 200 — 800 Myr. 

Of course, any application of our results to the Milky 
Way will require a more detailed understanding of how ver¬ 
tical modes couple to modes in the disk plane. As a satel¬ 
lite galaxy passes throug h the disk, it excite s both bend¬ 
ing and breathing modes dWidrow et al.|[2014ll . The former 
have been implic ated in the generation of galactic warps. 
iDebattist^ (l2014l l has shown that spiral arms generate com¬ 
pression and rarefaction in the disk but the reverse seems 
plausible, namely that the compression and rarefaction per¬ 
turbations due to a satellite, sheared b y differential rotation 
of the disk, generate spiral structure. Sellwood fc Carlberd 
ll2014h have argued that spiral activity in disk galaxies arise 
from the superposition of transient unstable spiral modes. 
Satellite galaxies and dark matter subhalos would seem to 
provide a natural seed for these instabilities. 


APPENDIX A: MATRIX ELEMENTS FOR THE 
FLUID CASE 

In this appendix, we describe the calculation of the matrix 
elements Mjk in Eg. 1161 As indicated in the text, we multiply 
both sides of Eq.[T0]by —ijjj and integrate with respect to 2 
from —00 to 00 . The left hand side becomes G^Cj. 

The first two terms on the right-hand side of Eg. llOl can 
be combined to give 
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( dpi \ 

y dz'^ dz dz J 


(Al) 


= J^^du (A2) 

= £ du (1 - n^) ((1 - n^) P,) (A3) 


where Afjk = k{k + 1) NjNk- We can then use Legendre’s 
differential equation as well as the identity 

(1 - u^) P'iu) = j iPj-iin) - uPO (A4) 

to write this expression in terms of the integrals 

= f du Pj{u)Pk(u) , (A5) 


= j duu^Pj{u)Pk{u) (A6) 

and 


In order to show that this is indeed a true mode of the 
system, we calculate the density directly from the distribu¬ 
tion function: pi = f dv fi. As discussed in the text, this 
integral formally diverges but can be handled by consider¬ 
ing a family of equilibrium distribution functions f\{E) that 
approximate fo but are continuous at if = 1. We then have 

J dvv^l^ = J dvv^ = - J dv fx = -po (B4) 

where at the last equality, we let A —>■ oo. Moreover, the 
term proportional to z^ vanishes. To see this, write dfo/dE — 
{l/2z) dfo/dz and perform the integration over velocities be¬ 
fore differentiating. The net result is that 


3iVi 1 
Stt 2-tj2/4 ’ 


(B5) 


which gives w (lAntonovI Il97ll : iKalnaisI 1 19731 : 

iFridman et al.llT98ll . 

The next mode (j = 3) is 


= J duuPj{u)Pk-i{u) , (A7) 

which can be ev aluated analytically (lArfken fc Webeill2005l : 
IWeissteinll2014l '). 

The third and fifth terms of Eq.[TD] can be combined to 
yield 


V>3 = Na (PVz) - P2{z)) = ^ (5z" - 1) {z^ - 1) . (B6) 

O 

Writing out in terms of trig functions 

^ (5£;2cos^6i-6£;cos2 6 l-bl) . (B7) 

8 

With the help of various trigonometric identities, we find 



= -4E-^^^(4-4) (A9) 

k 


while the fourth term gives 

dz^j^^dz = 4j2kNjNk {1% - 1%) . (AlO) 

APPENDIX B: HOMOGENEOUS SLAB 

For the lowest order even parity mode {j = 1) the potential 
can be written as 


V-i = ^ (■Bcos^ 6 >-1) (Bl) 

while density is constant: pi = 3Ai/47r. Since we are consid¬ 
ering an even parity mode, we can write the Fourier series 
for i/; as 


1*3,4 = 


7A3 

8 8 


and 


7A3 /5E' 


-3E 


'i/’3,2 = 


Thus, for the distribution function, we have 


(B8) 


(B9) 




21 * 3,2 


dE \2-uj^/4: 


cos 29 -|- 


2'i/'3,4 


2 -a;2/16 


cos 46 


(BIO) 

A straightforward, but tedious calculation similar to the one 
performed above leads to an expression for p3, which, when 
combined with Eq. [B7l yields a quadratic equation for 
whose solutions are uj/Q.^ = 2.07, 3.83. 


APPENDIX C: DYNAMICS IN THE COMPLEX 
ENERGY PLANE 


00 

'/’i = I] 1*1,n cosnd . 

n=0 


(B2) 


For the Spitzer sheet, we can find the density as a func¬ 
tion of the potential by integrating the distribution function 
(Eq.|43]) over v. 


We note that the n = 0 term does not contribute to the 
distribution function. The relevant Fourier coefficient is then 
i/’i,2 = 3Y1E/4 distribution function is then 


_ dfo 3Al E cos 26 
^ dE 2 2-a;2/4 ’ 


(B3) 


/ OO 1 

dvfiE) = —e-P. (Cl) 

The Poisson equation d'^ip/dz'^ = 2exp{—ijj) can then be 
integrated to obtain as expression for the force as a function 
of the potential: 
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The period for a particle with energy E is 
dv' 


1/2 


T{E) = -4 f 

'' 4^ma> 

^23/2^1/2 f 
Jo 


FW 

/2 


COS ifidip 


_ ^ — E cos^ 


1/2 


(C2) 

(C3) 

(C4) 


where i> max is the m aximum velocity of the particle along 
its orbit JArakilllflSRl l. In deriving the second expression, we 
write V and ip in terms of the parameter ip: 


V = {2E)^^^ sirup ip = E cos^ p (C5) 

More generally, we can write t along the first quarter of 
the particle’s orbit as a function of p: 


t {p, E) = 
Similarly 


£;\1/2 cos p'dp' 


ri 


,E) = Ef 
Jo 


_ ^ — E cos2 ^ 

sin ip'COS dip' 


For the lowered Spitzer sheet, p{ip) is given by 


(C6) 


(C7) 


pW) = fwe ^ (^(27r)^/^ erf (y) - 2®/^ye (C8) 

where y = {W — ipY^^ /n. As before, we integrate the Pois¬ 
son equation and find 


where yo 


FiiP) = { 87 r)F^ {Eiyo)-F{y))F^ 
and 


(C9) 


T{z) 


(27r)^/^e’""erf«-2®^^ 



(CIO) 


Finally, we write the complex Fourier coefficients of the 
basis functions ipj as 


1 df 

iPr^AE) = ^ dp — ip, {z [E, e)) (Clf) 

where z, t, and 6 are complex functions of E and of the real 
integration variable p through Eas. lCSI and EH 
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